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Abstract 



We study the effect of a Jovian planet on the gas distribution of a protoplanetary disc, 
using a new numerical scheme that allows us to take into consideration the global evolution 
of the disc, down to an arbitrarily small inner physical radius. We find that Jovian planets 
do not open cavities in the inner part of the disc (i.e. interior to their orbits) unless (a) 
the inner physical edge of the disc is close to the planet's location or (b) the planet is 
much more massive than the disc. In all other cases the planet simply opens a gap in 
the gas density distribution, whose global profile is essentially unchanged relative to the 
one that it would have if the planet were absent. We recognize, though, that the dust 
distribution can be significantly different from the gas distribution and that dust cavities 
might be opened in some situations, even if the gas is still present in the inner part of the 
^ ; disc. 

Concerning the migration of the planet, we find that classical Type-II migration (with 
speed proportional to the viscosity of the disc) occurs only if the gap opened by the planet 
is deep and clean. If there is still a significant amount of gas in the gap, the migration 
of the planet is generally slower than the theoretical Type-II migration rate. In some 
situations, migration can be stopped or even reversed. We develop a simple model that 
reproduces satisfactorily the migration rate observed in the simulations, for a wide range 
of disc viscosities and planet masses and locations relative to the inner disc edge. Our 
results are relevant for extra-solar planetary systems, as they explain (a) why some hot 
Jupiters did not migrate all the way down to their parent stars and (b) why the outermost 
of a pair of resonant planets is typically the most massive one. 

Key words: (stars:) planetary systems: formation - (stars:) planetary systems: pro- 
toplanetary discs - accretion, accretion discs - Solar System: formation. 
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1 Introduction 



The planet-disc inter actions have been subje ct of an increased interest since the discovery 
of the first exoplanet (IMayor fc Quelozlll995l ). Indeed, the first extra solar planets discov- 
ered were giant gaseous planets orbiting surprisingly close to their parent star, so that they 
are called h ot Jupiters. Accord ing; to classical planetary formation models in protoplan- 
etary discs (jPollack et al.lll996l ). hot Jupiters could not form where they currently orbit, 
because there was not enough solids to build a massive core that would accrete a gaseous 
atmosphere (there is much more solid material beyond the so-called 'snow-line' where wa- 
ter condensates into ice). In addition, gas heated by the star is more difficult to capture 
for the solid core than cold gas. It has been shown by lBodenheimer. Hubickyj. fc Lissauer 
( 120001 ) that in situ formation of hot Jupiter is not impossible, but it is not the most likely 
scenario. This suggests that the planets migrate in the disc, via angular momentum 
exchanges with the gas. 

In fact, the planet always exerts a positive torque on the part of the disc outside its or- 
bit, and a negative one on the part o f the disc inside its orbit (IGoldreich fc Tremaindll979 ; 



Lin fc Papaloizoulll979l ; IWardlll997l ). Reciprocally the disc exerts the opposite torques on 
the planet ; if the disc density profile is not perturbed by the planet, the sum of the torques 
is not zero but results in the so-called differe ntial Lindbla d Torque, which is negative and 
responsible for the inward type I migration (I Ward! 1 199 71 ). However, the torques exerted 
by the planet on the disc tend to repel the gas away from the planetary orbit ; whether 
the internal stress in the disc is strong enough to counterbalance this planetary torque 
and spread the gas into the void regions determines whether the corotation region of the 
planet is depleted or not, leading to the openin g of a more or less clean gap . For more 
detail on the gap opening process, we refer to ICrida. Morbidelli. fc Massetl (120061 ) and 
references therein ; the basic idea is that the density profile evolves in response to the 
presence of the planet, in order to achieve an equilibrium among three torques exerted on 
a fluid element. These torques are due to the gravity of the planet and to the viscosity 
and the pressure in the gas ; the two latter ones depend strongly on the gas density profile. 
The stronger is the planetary torque (it is proportional to the planet mass), the steeper 
has to be the density profile, and thus the deeper is the gap. 

If a clean gap is opened, the planet is locked into it because the outer and the inner disc 
are repellent for the planet. Thus, the planet is supposed to follow the viscous evolution of 



1974, hereafter LP74) ; this is the definition of type II migration. However, iQuillen et al 



Lvnden-Bell & Pringle 


However, 


Quillen et al. 



( 120041 ) noticed that if the planet is more massive than the inner disc, the outer disc is not 
able to push the planet inward at the speed of the natural gaseous accretion on to the star. 
Indeed, in that case the planet has more angular momentum to lose than the unperturbed 
disc would have, so that it migrates more slowl y. Meanwhile, nothing prevents the inner 
disc from accreting on to the star. Even more, IVarniere et al.l (120061 ) suggested that, as 
the planet exerts a negative torque on the inner disc, it accelerates the accretion of the 
latter on to the primary. Consequently, the inner disc accretes on to the central star faster 
than it would naturally. This makes the inner disc disappear before the planet migrates 
inward. The result is the opening of a cavity : a great depletion of the disc from the star 
to the planet's orbit. 
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Once a cavity is opened, the outer disc still tends to accrete on to the central star 
and to push the planet inward, while the planet can not take angular momentum from 
the inner disc anymore. A big mass of the planet with respect to the disc may slow down 
this process, but the final result should ineluctably be the inward type II migration of the 
planet toward the star. Reality, however, is not so simple. In fact, in this work we will 
show that, under some conditions on the gas disc viscosity, the orbit of the planet is not 
'gas-proof and viscous accretion can happen without forcing the migration of the planet. 

The cavity opening process and its feedback on planet migration are not easy to study 
and quantify with the help of numerical simulations. Indeed, th e computation of the inner 
disc evolution is prohibitive with classical hydro-codes. Thus, iQuillen et al.l (120041 ) only 
made the simulation of a planet o rbiting; inside a pre-ex isting cavity - showing that it did 
not migrate inward very fast - and lVarniere et al.l (120061 ) performed only a few simulations 
showing the opening of a cavity by a Jupiter-mass planet in a low density disc. The cavity 
opening appeared to be more rapid than the viscous time-scale, confirmin g the idea that 



the pl anet 'helps' the accretion of the inner disc on to the star. However, IVarniere et al. 



( 120061 ) did not provide an exploration of the parameter space, and could not study the 
planetary migration because the planet was about 100 times more massive than the disc. 
With the help of a new numerical scheme that simulate s efficiently and realistically the 
inner disc evolution ICrida. Morbidelli. &: Massetl (120071 ). we can perform such a study 
for the first time. This is the main scope of this paper. This analysis is particularly 
interesting because observations suggest that some protoplanetary discs most likely host 
cavities. Understanding what kind of cavity a planet can open in terms of size, shape and 
depth, in which conditions and with which life-time, may provide some keys to interpret 
these observations. We also stress that the ALMA interferometer will enable us to make 
precise images of the inner part of protoplanetary discs. In this frame, it is relevant to 
study the features resulting from planet-disc interactions that might be observed. 

In Section (2J we will shortly review the observations of protoplanetary discs that 
seem to host cavities. After a brief description of the algorithm used for our numerical 
simulations (Section EJ), we present the results in Section H] in terms of cavity opening and 
migration. These numerical experiments are interpreted in Sections [5] and [6] with a semi- 
analytical model. We discuss in particular the survivability of a cavity opened by a planet 
with regard to type II migration, and we find that the corotation torque can prevent this 
inward migration. In Section [8] we discuss the difference between gas cavities and dust 
cavities. Finally, we conclude and discuss the implications of this work in Section [9j 



2 Evidences for cavities in protoplanetary discs. 

Discs of gas and dust around young stars are heated by the star and by internal viscous 
dissipation. Consequently, the dust becomes an infra-red source. In the spectral energy 
distribution (SED) of the star, an infra-red excess reveals the presence of such a disc. 
The near-infra-red (NIR) emission corresponds to the hottest part of the disc, which is 
the closest to the central star, while the mid-infra-red emission is due to colder dust, 
corresponding to the outer part of the disc. Consequently, a lack of infra-red excess in the 
NIR part of the spectrum can be interpreted as the presence of an inner cavity in the disc 
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(jBeckwithlll999l ). Such a cavity discovered by SED analysis is sometimes called a spectral 
hole. 

One of the most clear examples of spectral hole is given in fig. 2 in lForrest et al.l (120041 ) 
for the T-Tauri star CoKu Tau/4 observed with the Spitzer telescope. The size of the 
spectral hole is such that there seems to be no dust at temperature larger than 123K. The 
authors convincingly deduce from this SED that dust grains in the disc of CoKu Tau/4 are 
excluded from the innermost 10 au ; this makes a remarkable dust cavity. The question 
whether this c avity could be opened and maintained by a planet has been addressed in 



Quillen et al.l (12004] ) . They suppose that a planet opened a gap in the disc at 10 au, and 
that the inner disc subsequently accreted on to the star. Then the age of the system 
has to be bigger than the viscous time-scale of accretion on to the star of the inner disc. 
This gives a lower limit for the viscosity. Given the viscosity, a minimum planet mass is 
required to open a gap. For the planet not to migrate immediately after the depletion 
of the inner disc, the planet has to be more massive than the outer disc, so that the 
angular momentum taken by the latter from the planet is a negligible frac tion of the 
angular momentum of the planet. From these considerations, iQuillen et al.l performed 
a simulation with a planet of mass 3 x 10 _4 M* on a circular orbit at 10 au in a disc 
with Reynolds number 10 5 and density S(r) = 10 _6 (l.lr p /r) M*au~ 2 for r > l.lr p and 
100 times smaller for smaller radii. The planet efficiently maintained the cavity for the 
duration of the simulation (100 orbits). However, they did not make a simulation of the 
process of cavity opening. 

The T-Tauri stars TW Hya, DM Tau and GM Aur also present a spectral hole. To 
interpret their SEDs, ICalvet et al.l (120021 . 120051 ) used the following model. They divided 
the disc into 3 components : an optically thick outer disc, its inner edge (represented as a 
wall directly exposed to stellar radiations), and an optically thin inner region. Adjusting 
the free parameters (the wall radius and height, the outer disc scale height, its mass 
and viscosity) they obtained a satisfactory fit of the SED of these three objects. They 
concluded that TW Hya and DM Tau present a 3-4 au wide cavity, while GM Aur should 
have the tru ncation radius of th e outer disc located at 24 au from the star. Previous work 
on GM Aur (IBergin et al.ll2004j ) based on different data, concluded that the cavity is 6-au 
wide. Such a ca vity could be ma intained by a 1.7 Jupiter mass planet on a fixed circular 
or bit at 2.5 au dRice et aT1l2003h. 



Dullemond. Dominik. fc Nattal (120011 ) proposed a different model for the SED of Her- 



big Ae stars that does not involve the presence of a planet. In their model, the inner part 
of the disc is removed, so that the inner rim is directly exposed to stellar radiations. It 
is thus heated a lot, and gets puffed up. This has two consequences: first, this rim is a 
strong IR source, leading to a peak at 2-3 \im in the SED ; second, the disc behind the 
puffed up inner rim is in the shadow. Thus, this region remains cold and does nearly 
not radiate in IR. If the disc is flared, however, the outermost part of the disc can be 
lightened by the star, and then becomes an IR source. In that case, the SED presents a 
gap between the peak at 2-3 fim due to the hot rim and the emission of the outer disc 
at larger wavelength, because it is colder. This is in good agreement with the SED of the 
Herbig Ae star AB Aur, with a rim of height 0.1 au at 0.52 au (which corresponds to the 
dust sublimation radius). 

More recently, interferometric observations at Plateau de Bure Interferometer clearly 
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show ed a remarkable dust depletion in the inner 50 au of the disc of LkCal5 (IPietu et al. 



20061 ) . The possibility that a massive planet or a brown dwarf is responsible for this is 



work in progress and will be addressed in a forthcoming article. 

Thus, cavities are relatively common features in protoplanetary discs. It is not easy 
to constrain their characteristics, as sometimes different dust distributions can give indis- 
cernible SEDs, but the evidence for cavities in protoplanetary discs is quite strong. 



3 Numerical simulations setup. 

Classical hydro-codes for the simulation of planet-disc interactions represent the disc with 
a 2D polar grid extending between an inner and an outer radius (Rf^ and -Rout)- The grid 
covers the planetary region and, provided that the resolution is good enough, satisfactorily 
describes the local planet-disc interactions. The grid cannot reasonably cover the full disc, 
for the following reasons. First, a prohibitive number of cells would be required to extend 
the grid from the planet region (a few astronomical units from the star) to the physical 
outer edge of the disc (often located at several hundreds of au from the star). Second, 
and more important, the simulation of the inner part of the disc requires a shortage of 
the time-step and an increase of CPU time roughly proportional to fi m a X Rfn ^ i n 
most codes, where Q is the angular velocity. Because the grid covers only a portion of the 
disc, the global evolution of the disc cannot be realistically described. Yet, it is the global 
evolution of the disc that governs type II migration. Thus, the accuracy of the type II 
migration observed in simulations can be questioned. Another point related to the global 
evolution of the disc is the evolution of its innermost part and its accretion on to the star, 
that is the opening of a cavity. Consequently, it seems that classical hydro-codes are not 
adapted for the study of these two problems on which we wish to focus. We developed 
an improvement to the classical algorithm, that we briefly present below. 



3.1 Code description 

The code used for this work is based on the hydro- code FARGO flMassetl l200o3 lbl) . To 
describe correctly the evolution of the disc outside the 2D grid, we add a ID grid extended 
all over the physical disc, from its inner edge (the radius of the central star or the X-wind 
truncation radius) to its outer edge. This non-azimuthally resolved grid is coupled to the 
2D grid at the boundaries of the latter in a conservative way. Thanks to this coupling, 
the disc perturbations due to the planet - computed locally in the 2D grid in FARGO - 
influence the global evolution described in the ID grid ; this evolution in turn affects the 
computation in the 2D grid as it provides realis tic, time-evo l ving boundary conditions for 
the latter. This numerical scheme, described in lCrida et al.l ( 120071 ). enables us to simulate 
the type II migration with an excellent reliability as well as the accretion of the inner disc 
on to the star, for a negligible additional computing cost with respect to a classical 2D 
hydro- code. 
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3.2 Units and parameters 



The central star mass is set equal to M Q , which is our mass unit. As length unit, we adopt 
the astronomical unit (au), and we set for simplicity the gravitational constant G = 1 (so 
that one year lasts 2tt time units). 

The initial density profile corresponds to a disc that evolved for some time under 
the effect of its own viscosity and was provided by iGuillot fc Huesol ( 120061 ). It can be 
approximated by : 

E(r) « 10" 5 exp(-r 2 /1320). (1) 

This initial density corresponds to the Minimal Mass Solar Nebula (jHayashil Il98ll ) at 
the location of Jupiter (5.2 au). The aspect ratio H/r is constant in space and time; its 
default value is 0.05. These parameters (density and aspect ratio) will be changed in some 
test runs to study their influence on the results. The sound speed is c s = HQ, where H 
is the height of the disc. The equation of state is isothermal (the pressure is P = pc 2 , 
where p is the volume density, so that it becomes P = Y>c 2 in our 2D formalism, after 
integration on the disc's scale height). 

If not specified otherwise, in most of the runs the 2D grid extends radially from 
Rf® = 1.75 to R^ t = 15. au and is divided in N r = 165 elementary rings, themselves 
divided into N s = 320 sectors. The ID grid extends radially from 0.58333 to 100 au, over 
iV r 1D = 1193 elementary rings. 



4 Results 

We computed a few simulations with a Jupiter mass planet (q = M p /M* = 10~ 3 ) initially 
placed on a circular orbit at r p = 5 au in a disc with different viscosities. The kinematic 
viscosity of the gas, v, is constant in space and time. Its values in the simulations are 
such that the Reynolds number at the initial location of the planet (1Z = r 2 Ql v jv) goes 
from 10 3 ' 8 to 10 6 . 



4.1 Density profile and cavity opening 

Fig. [1] shows with bold lines the density profiles for four simulations after half a viscous 
time (t u = ur~ 2 ), or only a quarter of viscous time in the least viscous case. The depth 
of the gap opened by the planet strongly depends on the Reynolds number. This is not 
surprising, as it i s well-known that viscosity plays against gap opening (see for instance 



Crida et al.l 120061 . and references therein). The gap is centred on the planetary orbit, 
which is not any more at 5 au because the planet has migrated. 

Concerning the opening of a cavity, we see that the density in the inner disc is signif- 
icantly smaller than in the outer disc. In particular the least viscous cases show an inner 
disc with very low density, so that a sort of cavity is formed : in fact, the density is about 
5 times smaller for r < 2.5 than for r > 3. The profile clearly shows a wall at r pa 2.5, 
beyond which the disc profile is about flat while the density is negligible at the base of 
the wall. 

The thin dotted lines show the density profile obtained in similar simulations, at the 
same time, in the same conditions, but without the planet. The density also decreases in 
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Figure 1: Bold lines: Surface density profile at t u /2, for different Reynolds number, 
labelled on the corresponding panel. For TZ = 10 55 , the time is t u /4. Thin dashed lines: 
surface density profile at the same times for a similar disc with no planet. The units of 
mass and length are those specified in Section 13721 : solar mass and astronomical unit. 
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Figure 2: Migration of a Jupiter mass planet for different Reynolds numbers of the disc. 



the innermost region of the disc. So, it seems that the observed depletion of the inner 
disc is not fully caused by the planet, as will be discussed in Section [51 



4.2 Migration rate 

Fig. [2] shows the evolution of the semi major axis of the planet with time, in units of 
the viscous time t u . In proper type II migration, the migration rate is proportional to 
the viscosity and thus, in this time unit, it should be independent of 1Z. This is indeed 
the case for TZ ^ 10 5 : the three curves are almost linear and overlap, at least in the 
first part of the evolution. However, for higher viscosities, a very different behaviour is 
observed. When TZ decreases below 10 5 , the migration rate becomes slower than expected. 
As viscosity increases, the planets migrates inward more and more slowly with respect to 
the viscous time, and the migration is even stopped for TZ = 10 415 . We will refer to this 
case hereafter as the stationary case. For larger viscosities, the migration is reversed and 
the planet moves outward. 

This result is particularly surprising and new. This outward migration is not an effect 
of the resolution of the grid; we have recomputed the simulation with TZ = 10 4 , using 
N r = 330, N s = 640 : the planet also migrates outward, at about the same speed (see top 
panel in Fig. [3]). With this high resolution, the size of a cell around the planet is about 
a tenth of the Hill Radius in radius and a seventh of the Hill radius in azimuth, so that 



8 



5.8 
5.6 
5.4 
5.2 
5 
4.8 
4.6 

5.4 
5.3 
5.2 
5.1 
5 
4.9 
4.8 
4.7 
4.6 



I 


I I 




' R = 10 4 




N r = 330 ; N s = 640 


I 


N r = 165 ; N s = 320 



0.5 



0.5 



1.5 



R = 10 415 
H/r = 0.04 

H/r = 0.05 

H/r = 0.06 



1 

time [t v ] 



1.5 



Figure 3: Top panel : Outward migration of a Jupiter mass planet for 1Z = 10 4 and two 
different resolutions. Bottom panel : Quasi stationary evolution of a Jupiter mass planet 
for 1Z = 10 415 and 3 different aspect ratios. 



the Hill sphere of the planet is covered by about 220 cells. This is largely sufficient; in 
particular the corotation zone of the planet, which plays a crucial role as will be shown 
further, is correctly simulated. We also investigated the effect of the aspect ratio on the 
migration rate, from the stationary case. It appears that this parameter only plays a 
marginal role (see bottom panel in Fig. [3]). 

One may wonder if the stationary case observed for TZ = 10 for a Jupiter mass 
planet is a feature valid for any planet mass. Fig. H] shows that it is not the case. The 
more massive is the planet, the faster it migrates inward, approaching the classical type II 
regime. For planets lighter than Jupiter, the migration is directed outward. However, we 
will show in Section [6] that if the viscosity is independent of radius, there is a stationary 
orbit further from the star which these low mass planets tend to, asymptotically. 

5 Modelling the density profile 

In this Section we build a model to understand the shaping of the disc density profiles 
observed in Fig. HJ Our approach is done in three steps. First, we review the global 
viscous evolution of a disc that is not perturbed by any planet; then, we discuss how a 
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Figure 4: Migration of giant planets of different masses in a disc with K = 10 415 and 
H/r = 0.05. 



gap is opened by a massive planet and finally we analyse in which situations the resulting 
density profile may show a significant depletion (a 'cavity') in the full inner part of the 
disc, up to the planetary orbit. 



5.1 Global disc's viscous evolution 

The natural viscous evolution of a gaseous disc is described in LP74. Using their partial 
differential equation (14), one can compute, with a step by step integration, the evolution 
of any initial profile. If the initial profile has a Gaussian shape (£(V) = Eoexp(— ar 2 ), 
where a is an arbitrary constant), there is an explicit solution, given by equations (18) 
and (18') of LP74. The initial shape is conserved, the surface density being: 



E(r) = E T- 5/4 exp 



ar 



T 



where T = Ylavt + 1. 

This solution is valid only for a disc extending between and infinity in radius. How- 
ever, the inner radius of the disc, -Rmf , is never ; it is at least the radius of the central 
star, most likely the corotation radius at a few tenth of au where the lines of the magnetic 
field reconnect (the so-called X-point; Shu et al. ( 19971 )). It could even be bigger, in the 
case o f Jet Emitting Disc (J E D) : indeed, a t the b ase of the jet, which may be sever al-au- 
large feacciotti et al.l fbOQSh : ICoffev et all (120041 ) ). accretion is dominated by the torque 
exerted by the jet and is much larger than the standard accretion ; so, the outer radius 
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Figure 5: Surface density profile of an unperturbed disc at t — t u /2. Plain lines: profiles 
from simulations with TZ = 10 4 , 10 4 ' 5 , and 10 5 . Dashed line: profile from equation ([2]). 



of the jet could be considered as the inner open boundary of the standard disc. The 
explicit solution of the LP74 equations with a finite inner edge of the disc comes from 
equation (25) of LP74. With our units (G = M* = 1), it gives: 

WM) = EoT"^ (^^) exp (-£) • (2) 

Thus the density starts from at r = i?i n f and grows with r until r ^> R- m t, where a 
classical Gaussian shape is reached. This kind of profile is illustrated by the thin dashed 
lines in Fig. [IJ Fig. [5] shows the agreement between the disc profiles obtained for various 
viscosities in numerical simulations and the theoretical profile given by equation (j2J) with 
S = 1.224 x 10~ 5 and a = 1/1320 at t v /2, (id est T = 1.114). All profiles overlap with 
each other. This shows that, once the time is renormalised relative to the viscous time, 
the evolution of the disc density distribution is independent of the Reynolds number. 
Also, it shows the excellent agreement between equation ([2]) and the numerical solution 
achieved with our simulation scheme. The little discrepancy visible in Fig. [5] comes from 
the fact that the original profile in the simulations is not given by equation (j2J) at T = 1 
but by equation (JTJ). 



5.2 Gap opening in an evolving disc 



Equation (14) of ICrida et al.l (120061 ) provides a way of computing semi-analytically the 
density profile of a gap for any disc and planet parameters. This equation comes from 
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the equilibrium between the torques due to the gravity of the planet, the viscosity and 
the pressure in the disc. It applies outside the corotation zone of the planet, which 
approximately extends on each side of the orbit over a width xf ~ x~ ~ 2Rh, where Rh 
is the Hill radius of the planet, the equation reads: 



R H dE 
E dr 



3 fl 

0.35g 2 ftr p 3 rA- 4 sgn(A) - -— 



H 



r 3 
— Ka" + 



n 



(3) 



r r . 



IRh Vt 



where a" 



Rh 

Notice that a 



-1.2 



+ 200 



Ri 



-10 



and A 



r — r. 



p- 



boundary condition is needed to solve this differential equation with a 
step by step procedure. The solution of this equation corresponds to the gap profile E(r) 
in an ideal disc in steady state, under the influence of a non-migrating planet. Without 
the presence of the planet, the unperturbed profile of such a disc would be proportional 
to r -1 / 2 because v is assumed to be constant over r. Thus, to compute the profile in 

Eo/V-Rmf & t r = -^inf an d we compute the 
For the outer disc, we start with E 



the inner part of the disc, we start from E 



density step by step until r~ 



x 



R* 



at r = -R sup and follow equation 



J 0/ Y -'•-sup 

I) down to r + = r p + xf. Unfortunately, T^(r~)/y/r~ 
is not necessarily equal to E(r + )/-y/r+. Thus, we choose the lowest of these two values, 
say the one at r~, and adjust the outer edge of the gap at r + + 5, with 6 such that 
E(r + + 5)/\/r + + 5 = E(r~)/Vv ■ As shown in ICrida et al.l ( 20061 ). this procedure gives 
a satisfactory approximation of the gap depth. 

The use of equation ([3]) as explained above provides the gap profile in a disc whose 
unperturbed profile is E(r) = E r -1 / 2 . Thus, multiplying this gap profile by y/r/T^o gives 
the gap profile in terms of fraction of the unperturbed profile. Let us denote cr(r) this 
fractional profile. We think that it is reasonable to assume that the fractional profile is 
independent of the unperturbed profile. This is supported by Fig. [TJ which shows that 
the planet opens a gap in the natural profile of the disc, without essentially changing it. 
Consequently, as a simple model of the gap profile in an evolving disc, we suggest that 
the unperturbed profile given by equation (T5]) can be multiplied by the fractional profile, 
obtaining : 

E(r) = a(r)E LP74 (r,t) . (4) 

Fig. O displays this model in the 1Z = 10 415 case. The secondary panel shows the 
fraction profile cr(r). In the main plot, the dashed line shows the unperturbed profile 
Elp74(?", t) at t = t u /2. The gap profiles from the numerical simulation at the same 
time (thin line) and from equation (jlj) (bold line) are compared. The agreement is quite 
satisfactory, in particular if one keeps in m ind that we o nly h ave a semi-analytic estimate 
of the gap profile. Indeed, as discussed in ICrida et al. fl2006h . equation ([3]) provides gap 
profiles that are a bit too narrow and shallow for TZ ^ 10 5 (while it provides slightly 
too wide and deep gaps for 1Z > 10 5,5 ). This is observed here as well. However, the 
gap edges in the model correspond quite well to the ones from the simulation in terms 
of position and slope. In addition, the profiles in the inner and outer disc nearly match 
those observed in the simulation. 
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Figure 6: Gap profile in an evolving disc with 1Z = 10 415 : the thin line results from 
numerical simulation while the bold line comes from our simple model. Dashed line : the 
unperturbed profile £lp74(?~, t„/2) from equation (j2|). Secondary panel: 'fraction profile' 
cr(r) of the gap deduced from equation 
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We present here only this case, corresponding to a non- migrating planet. In the cases of 
a migrating planet, the situation is similar. This shows that this simple idea of multiplying 
the unperturbed profile by a schematic gap profile is valid in first approximation. As a 
matter of fact, a slight depletion of the inner disc with respect to the unperturbed profile 
is visible for cases with inward migrating planets. We will explain this in Section [7J after 
having developed a model that reproduces the migration rates. 



5.3 Cavity opening 

As we have seen above the disc density profile can be described effectively as the product 
of the fractional gap profile with the profile described by equation (j2J). This equation 
gives profiles which are proportional to (l — (r/.R m f) _1 ' 2 ) exp(— ar 2 /T). The first term is 
about 1 for r > 30.R m f, so that the profile at large radius is about Gaussian. However, 
for r/R [n f < 5, the first term shapes the profile growing from at r = _R in f , with a slope 
proportional to (r /R in f)~ 3 / 2 . Thus, it seems that R\„f is a key parameter fo r the profile of 



the innermost part of an evolving disc (see also iLubow &: D'Angeloll2006l ). 

Fig. [7] shows the density profiles at t — t u /10 for discs with 1Z = 10 5 and four different 
values of -R m f, as they result from numerical simulations. The plain line corresponds to 
-Rmf = 0.65 au as before (nominal case). The dotted line corresponds to a four times 
smaller value of -R m f. All other parameters are the same in the two simulations. In 
both cases, the Jupiter mass planet has opened a gap in the disc and migrated down 
to r v m 3.85. The outer disc profiles almost overlap. The inner disc profiles, however, 
present a huge difference : the maximum density is nearly twice bigger in the small R{ n f 
case than in the nominal i?i n f case. In the nominal case, the inner disc is being depleted 
and a cavity is appearing. In the small i? in f case, there is only a gap. The tendency to the 
opening of deeper cavities with larger i? in f is confirmed in the last two simulations (short 
dashed and long dashed curves, respectively, in Fig. [7]). 

In conclusion, depending on the width and shape of the gap and on the position of the 
planet r p with respect to the inner edge of the disc -R m f, the density profile may show a 
cavity or not. Indeed, if the inner edge of the gap falls in the flat profile zone (r > 10 -Rmf), 
a classical gap is shaped (case 1Z = 10 4 in Fig. [Hand case i?i n f = 0.16 in Fig. [7J). If it falls 
between and ~ 4 -Rmf, the density in the inner disc cannot grow enough before that the 
inner edge of the gap is reached and the planet seems to open a cavity (case 1Z = 10 5 ' 5 in 
Fig. [I]). Thus, the main parameter to determine whether a planet opens a cavity or not is 
the ratio between r p and -R m f. This should be taken into account in the interpretation of 
numerical simulations. To perform realistic simulations, it is necessary to have a realistic 



value of -R in f. Thus, our code with coupled 2D and ID grids (ICrida et al.l 120071 ) seems to 
be a tool of choice as it can handle arbi trarily small -R 



nf- 



Our result differs from the one by I Varniere et al.l (120061 ) . They claimed that the 



depletion of the inner disc was faster than viscous because of the negative torque exerted 
by the planet. However, a planet in equilibrium in the middle of the gap simply transfers 
to the inner disc the torque that it feels from the outer disc ; otherwise, the planet would 
move with respect to the gap. More precisely, in type II migration, the planet generally 
feels a total negative torque, and migrates inward, together with the disc and the gap ; 
thus, the negative torque that the planet exerts on the inner disc is a bit smaller in 
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magnitude than the one it feels from the outer disc. The difference corresponds to the 
loss of angular momentum of the planet. In the absence of the planet, there would be 
gas in the gap ; this gas would also feel a negative torque from the outer disc, exert a 
smaller negative torque on the inner disc, and migrate inward losing angular momentum. 
If the planet has the same mass as the gas that was initially in the gap, it has exactly the 
same effect on the disc. Consequently, we believe that the presence of the planet does not 
modify substantially the evolution of the inner disc - exce pt maybe if the plane t is much 
more massive than the disc, which was indeed the case in lVarniere et al. ( 2006 ) analysis. 
This is confirmed in our simulations (see Section [7]) . 



6 Modelling the migration rate 

It seems quite logical that, when the planet opens a clean gap, its migration follows a 
proper type II regime. With a Jupiter mass planet and a disc aspect ratio of 0.05, this 
happens for 1Z > 10 5 (see Figs [1] and [2]) . 

However, for smaller Reynolds numbers, the gap is not completely gas-proof. The gas 
in the gap has two major consequences : (i) it partially sustains the outer disc, effectively 
reducing the torque felt by the planet from the outer disc (ii) it exerts a corotation torque 
on to the planet. The possibility of gas flowing through the gap decouples the planet from 
the gas evolution. 

In this Section we show with a simple model that taking into account these effects 
allows us to explain the evolution of the planet as a function of the various parameters. 
Our model is based on previous works on the corotation torque flMassetllgOOlh the VISCOUS 
evolution of accretion discs (LP74), and the shape of gaps (jCrida et al.l 120061 ) . 



6.1 Classical type II torque 

In an accretion disc, the viscous stress is such that angular momentum flows outward 
while matter falls on to the star. In a Keplerian, circular disc with v and S independent 
of the radius, the torque exerted by the part of the disc extending from a given radius 
r to infinity on the inner part {r < r } is T v = — 37r£z/r 2 fi (it can be easily found 
from the strain tensor). It causes a mass flow of gas F, carrying the equivalent angular 
momentum: T v = Fr Q 2 Q Q = (27ir Ilv r )r 2 Q , where v r is the radial velocity of the gas. 
In this model, v r = — |^, which can also be found from the Navier-Stokes equations. 
This gives the following equality, which we will use further : 

" = — ^-r • (5) 

If a planet opens a deep gap in such a disc, no gas flow is allowed through the planetary 
orbit. The outer disc is maintained outside of the gap by the planet, and an equilibrium 
is reached so that the planetary torque balances T„. Consequently, the planet feels from 
the outer disc the torque T v . This torque is proportional to the viscosity and not to the 
planet mass. This is the case of standard type II migration. 

In a more realistic, viscously evolving disc, the scheme for type II migration is the 
same, but the above formula for T v is no longer valid. In that case, the equations of LP74 
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provide the density, the viscous torque, and the radial velocity as a function of radius and 
time. In our case of a disc with i? in f > 0, it gives : 



T v = 37n/£ T- 5 / 4 (/i-Mexp(-^ 



(6) 



£ L P74 = T v / (3-KVy/r) (7) 
BT 

F/ (27rrS LP74 ) . (9) 



F 



where h = r 2 Q = y/r is the specific angular momentum. Notice that equation (jBJ) is 
exactly equation (25) in LP74, while equation (J2j) is equivalent to equation (TSJ). 

Thus, in standard Type II migration, we consider that the planet feels from the disc 
a torque 

T n = Fh = 27rr + £ L p 74 (r + X(r + )v / r r , (10) 

where r + = r p + x s is the radius of the external edge of the gap, and £lp74 and v r come 
from equations (GO) and jHJ) respectively. 



6.2 Torque exerted on the outer disc by the gas in the gap 

The gas in the gap, the density of which is denoted S gap , exerts on the outer disc a 
positive viscous torque that is given by equation (flOl . with S gap instead of £lp74 and 
the opposite sign. This torque partially sustains the outer disc, and therefore needs to 
be subtracted from the torque that the planet would suffer from the outer disc if the gap 
were clean (given by equation ffTUl) ). So, denoting by / the ratio S gap /E L p74 we have 

T(i) = -/T n . (11) 

We now discuss how to evaluate / in practice. We have presented in Section [5] a way 
to compute semi-analytically the gap profile and the gap depth. However, making a step 
by step integration until the bottom of the gap it is not very convenient. Consequently, 
we looked for a simple empirical formula for the ga p depth as a f unctio n of the viscosity, 



the aspect ratio of the disc and the planet mass. ICrida et al.l (120061 ) showed that the 
density inside the gap is less than 10 per cent of the unperturbed value (i.e. / < 0.1) if 
and only if : 

7>=L!L + »<i. (12) 

Using equation we have computed the depth of the gap for various values of the 
parameter V. For each value of V, we impose q = 10~ 3 and H/r = 0.05, and find the 
corresponding viscosity. Then, we use these parameters in equation (jHJ) ; the obtained 
gap depth is shown as big dots in Fig. [HI We repeat the same operation for q ranging 
from 5 x 10~ 4 to 2 x 10~ 3 ; the results are reported as crosses in Fig. [HJ Furthermore, we 
impose q = 10~ 3 and v — 0, and find the corresponding H/r and the resulting gap depth. 
We repeat the operation for v ranging up to its maximum possible value compatible with 
V ; the results are shown as dots in Fig. [HJ Finally, for a given value of V, different 
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Figure 8: Gap depth (measured as ratio of the gap surface density to the unperturbed 
density at r = r p + 2R H ) as a function of V . The data points for each value of V are 
obtained from the integration of equation ([3]), assuming different values of v and H/r and 
keeping q = 10~ 3 (points) or different values of v and q and keeping H/r = 0.05 (crosses) ; 
the big dots correspond to the gap depths obtained for different values of v and keeping 
both H/r = 0.05 and q = 10~ 3 (see text for a more precise description of the sets of 
measures). The bold line is an approximate fit of the data. 



depths are observed because the shape of the gap is not the same whether the viscosity 
or the pressure dominates. However, there is a very clear tendency for deeper gaps with 
decreasing V, as could be expected. The bold line in Fig. [S] shows the fit that we adopt : 

J(^-°- 541 )/4 if V< 2.4646 
W j "\l-exp(-^) if V> 2.4646 ^ 

where f(V) stands for S g a P /SLP74- 

As the gap depth given by equation (j3J) is not very precise and shows variations at 
fixed V, this fit is clearly a crude approximation of the real depth. However, the evolution 
of the gap depth as a function of the parameters is satisfactorily described by f(V), in 
particular if one imposes q = 10~ 3 and H/r = 0.05, which is the most common case for 
us. As we are looking for a simple model, this approximation is sufficient for our purpose. 
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6.3 Corotation torque 



In iMassetl (120011 ). the corotation torque is computed as 

3 



T, 



c 



_ X a ^-Jno 3~ (Za 



where T,^ is the unperturbed surface density (assumed to be independent of radius), z s = 
(lZ/2ir) 1 ' 3 x s /r p , the function T{z) = z~ 3 — z~ i g(z) /g'(z) a nd q is a linear combination of 
the Airy functions Ai and Bi defined in equation (A18) o f IMassetl ( 20011 ) . This formula 
can easily be generalised for an arbitrary profile of £ (see IWardlll99ll ). and be rewritten 
as: 



T. 



c 



3 2 4 dlogg/g) ^ 

lip J/j £J 

4 d log r 



1 



where B is the second Oort constant (B 
4^((s/27t)<1/3)). a 

The function is shown in graphical form in Fig. 2 in 

(which is our case if we consider x s = 2Rh, q 
Thus, we may write : 

T 3 
J-c 



x s /r p ) 3 TZ / 
0/4 in Keplerian rotation), and T{z) 



Massetl (l200lh . For z < 0.1 

20z. 



20a;,fi„r„£z/ 
4 s P P 



10~ 3 and ft > 10 3 - 6 ) one has T{z) « 
dlog(£/B) 



dlogr 



In IMassetl ( 120011 ) . the approximation that z/ and S are independent on r is made. 
Thus, we can use equation (jSJ) and write equivalently : 



T c = — l0x K Q„r n Ev, 



v v 



dloggAg) 
dlogr 



(14) 



Here, £ is the density inside the gap, which is equal to SLP74(r + )/(7 : '), and v r is the radial 
velocity at r = r p , which is given in equation 0. The term dl ° g 1 ^ B - > is computed from 
the unperturbed density S L p 74 (r), assuming a Keplerian rotation of the disc. 



6.4 Total torque exerted on the planet 

The total torque felt by the planet is therefore , T p = Tu — Tu) +T C . From equations (fTUl) . 
( TTTj) and (fT4"|) . the total torque reads: 

T p = 27rr+S LP74 (r+)t; r (r+)v / r+x 

1 _ fCP) — 15 Xs Qprp2 Vr ( r p) f cp\ dl °s(^/ B ) ] l-^J 

J V ) 2nr+ ^+v r (r+) J y > dlogr 

with V and / defined respectively in equations (fT2"l) and (fT3"j) . while r + = r p +x s , x s = 2i?^, 
and Elp74 and v r come from equations (J6j) - (J9D - This expression involves (directly or via 
v r or SLP74) the viscosity of the disc u, its aspect ratio H/r, the planet to primary mass 
ratio q, the radius of the planetary orbit r p , and the radius of the inner edge of the disc 
i?i n f. In the following, we test it against numerical simulations for a wide range of these 
parameters. 
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6.4.1 Dependence on the Reynolds number 

In the top panel of Fig. [9] we plot the total torque felt by a Jupiter mass planet located 
at r p = 5 au in a disc with 0.05 aspect ratio and R in f = 0.65 au as a function of TZ. 
The torque given by equation (FIB]) with t = t u /100 is plotted as a bold line. The crosses 
with the error bars represent the torque felt by the planet in the numerical simulation, 
measured at t = t^/100 (so that the planet is still at about r p = 5 au). The error bars 
correspond to the maximum and minimum migration rates measured between t^/200 and 
tv/50. The torque Tn, corresponding to classical type II migration, is drawn as a thin 
line ; it is proportional to v oc TZ~ X . The bottom panel shows, as reference, the gap depth 
f(V) as a function of TZ. 

As one can see in Fig. [9j the model reproduces very well the general tendency. For 
TZ > 10 5 , f(V) < 0.125 and T p is close to T n , negative and proportional to the viscosity, 
as expected in proper type II migration ; in other words, Tn dominates T(j) and Tc in 
the calculation of T p . For lower Reynolds numbers, the effects of the gas in the gap 
counterbalance the classical T n . For TZ < 10 4 , f(V) > 0.7 and the total torque is positive. 
Both the measures from the simulations and the results of equation ffT5]) show a minimum 
of the torque felt by the planet at about TZ = 10 4 ' 5 . This shows that type II migration 
speed is bounded. One could think that an increase in the gas viscosity, and therefore 
of the gas accretion rate, would increase the type II migration speed of a giant planet. 
In reality, this favours the filling of the gap, which decouples the planet from the disc 
evolution. 

Notice however that for very low Reynolds number (high viscosity) our model reaches 
its validity limits. In fact, our model predicts a very fast outward migration, which is not 
observed in the simulations. In these cases, f(V) is more than 0.8 (essentially no gap), 
so that the basic idea on which our model is built is no longer valid. Eventually, for large 
enough viscosity, the planet has to migrate inward, at a Type I migration rate. 

6.4.2 Dependence on the planet mass 

In Fig. [101 the torque exerted on a planet placed at r p = 5 in a disc with TZ = 10 4 ' 15 , 
H/r = 0.05 and i?j n f = 0.65 is plotted as a function of the planet mass. The bold line 
corresponds to our model, and the crosses correspond to numerical experiments (with the 
error bars computed with the same prescription as before). Once again, we recover the 
observed trend. The discrepancy at very low mass comes in fact from the estimation of 
the gap depth f(P) ; for q < 5 x 10 -4 , the density in the 'gap' is larger than 0.8, so 
that the gap is shallow and / is not very accurate. Using directly equation ([3]) instead of 
prescription equation (fT3l) for the gap depth gives a better result, shown as a thin line; 
this validates again our model. 

The torque felt by the planet is a decreasing function of its mass. The explanation for 
this behaviour is quite easy: the lighter is the planet, the shallower is the gap, and the 
more important are the effects of the gas in the gap on the total torque. Symmetrically, 
the more massive is the planet, the deeper is the gap it opens, and thus the closer to 
type II migration is its migration regime. 

This result is relevant for extra-solar planets as it can explain why, of the many couples 
of resonant planets, typically the most massive object is the outermost one. Indeed, 
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tion [321 as a function of the Reynolds number. Bottom panel : the gap depth f(V) as a 
function of the same parameter. 
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resonant capture is possible only if the outer planet migrates inward faster than the inner 
planet which, according to Fig. [TUJ requires that the outer planet is the most massive. 

6.4.3 Dependence on planet accretion 

The planet mass is not necessarily constant. Planets are supposed to accrete gas contin- 
uously. Gas accretion by the planet perturbs the gas flow through th e planetary orbit ; in 



addit ion, mass accretion exerts an additional torque on the planet (ILubow fc D'Angelo 



20061 ). Therefore, in principle, accretion could ch ange t he dy namics. We have imple 



mented planetary accretion following the recipe of iKleyl (119991 ) : it consists in removing 
a fraction of the material in the Hill sphere of the planet and adding it to the mass of 
the planet. The accretion rate, expressed as a fraction of mass removed per time unit is 
imposed as an input parameter. More precisely, we apply the input removal rate / accr in 
the inner part of the Hill sphere (extended up to 0.45 Rh), and two thirds of / accr in the 
region from 0.45 Rh to 0.75 Rh- 

Considering 1Z = 10 3 8 and H/r = 0.05, we performed simulations of a Jupiter mass 
planet initially on a circular orbit at 5 au, with input removal rates / accr = 0.5Q P and 
/accr — 0.1Q p . The results in terms of migration and accretion are shown in Fig. [11] 
and compared with the no accretion case. Another simulation with input removal rate 
0.5Q p and initial mass of the planet 5 x 10 -4 is also performed. The results show that 
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accretion does not prevent outwards migration. The comparison between the top panel 
(migration) and the bottom panel (mass evolution) shows that the migration speed seems 
to be governed by the planet mass, and not by the accretion rate. While the accretion 
rate is constant with time in all cases, the migration rate is not, and its evolution follows 
the evolution of the planet mass. As expected from the discussion above, the outwards 
migration rate is slower for more massive planets ; if the planet reaches 3 Jupiter masses, 
its migration vanishes ; the migration is directed inwards for more massive planets. 



6.4.4 Dependence on the radii of the planetary orbit and of the disc inner 
edge. 

One may wonder about the dependence of the migration rate on the planet's location in 
the disc. This is governed by the radial dependence of viscosity and aspect ratio. 

If the viscosity is constant, the Reynolds number TZ increases as \fr. Thus, the depth 
of the gap opened by the planet increases with the planet's location r p . As a consequence 
a planet should migrate inward if it is sufficiently far from the star, and outward if it is 
sufficiently close. The migration paths therefore converge towards the stationary solution, 
at some specific radius, dependent on viscosity and planet mass. 

However, if the viscosi ty of the disc depends on radius as described in an a model 



( IShakura &: Sunyaevl Il973l ). the Reynolds number is independent of radius. If the disc 
is flared, H/r increases with r and therefore the depth of a gap opened by the planet 
decreases with r p . Consequently the planet tends to move outward if it is far from the 
star and inward if it is close. Any stationary solution would therefore be unstable. 

If the disc is not flared and the Reynolds number does not depend on r, then the 
behaviour of the planet depends on the ratio r p /_R inf . This effect is subtle. If the ratio 
r P /Rini is smaller, the planet is in a location of the disc where S and v r have a steeper 
positive slope, and v r is also larger in absolute value, following equation In our 
model, the gas radial velocity multiplies all the torques that appear in equation ( IT51) . 
Therefore, the only dependence on r p / i?i n f is through the gradient of the gas radial velocity 
v r (r p ) /v r {r + ) and the gradient of the logarithm of the density that appear in equation (fT5l) 
for the corotation torque. Fig. [12] shows as a bold line the torque felt by a Jupiter mass 
planet for 71 = 10 4 ' 15 at t — t^/lOO, as expected in our model. The gradient effect becomes 
important when Ri n {/r p ^1/2. The crosses with the error bars correspond to measures 
from numerical simulations. Once again, our model reproduces very well the observed 
trend. 

This result might be relevant to explain the existence of hot or warm Jupiters. Al- 
though several solutions have been proposed, the issue of why these planets did not 
migrate all the way into their parent stars remains open. Figure d2] shows that, for the 
parameters used in that calculation, a planet approaching the inner edge of the disk has 
eventually to stop at the location where the migration rate turns from negative to posi- 
tive. Therefore, the survival of some hot / warm Jupiters against migration seems to be 
a natural consequence of the local surface density gradient at the inner edge of the disk, 
predicted by LP74. Notice however that some hot Jupiters are so close to the central star 
that they are presumably well within the corotation radius. As the disk inner edge should 
not have been closer than the corotation radius, these planets somehow managed to mi- 
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grate past the inner edge of the disk. In the framework of our model this may be possible 
for massive planets in low viscosity disks, opening deep and clean gaps. In these cases, 
the planets presumably stopped migrating when their outer 1 : 2 mea n motion resonances 
reached the inner disk edge, as proposed in iKuchner fc Lecarl (120021 ) . In conclusion, the 
dynamical diversity of the exoplanets could be a consequence of the physical diversity of 
the protoplanetary discs. 



7 Effect of planet migration on the inner disc deple- 
tion 

In Fig. [6] we have shown that in the case of a non-migrating planet, the density profile in 
the inner part of the disc coincides with the unperturbed profile given by LP74 equations. 
This showed that the overall mass flow through the disc was unperturbed by the presence 
of the planet. 

Fig. [13] is the same as Fig. [6] but for a case where the planet moves inward at slow rate. 
In this case there appears to be a deficit in the inner part of the disc, while the outer part 
shows a little over-density, relative to the unperturbed case and to the prediction of our 
simple model. This reveals that only a fraction of the unperturbed mass flow effectively 
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Figure 13: Gap profile in an evolving disc with 1Z = 10 4 ' 5 : the thin line results from 
numerical simulation while the bold line comes from our simple model. Dashed line : the 
unperturbed profile £lp74(?", t u /2) from equation (|2D. Secondary pa nel: 'fraction profile' 
cr(r) of the gap integrated from equation (14) of ICrida et al.l (120061 ) . 
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passes from the outer to the inner disc. 

Based on the model developed above, we can illustrate this fact and develop an a 
priori estimate of the inner disc mass deficit. 

Due to the migration of the planet, the surface of the inner disc shrinks at a speed: 

A p = 27Tr p f p . 

In the unperturbed case, if we draw a fictitious boundary between the inner and the outer 
disc moving at the radial velocity of the gas, the surface of the inner disc shrinks at a 
speed : 

^unp 2,7TrpV r . 

Thus, for the inner mass of the disc to remain the same as in the unperturbed case, the 
mass flow into the inner disc through a planet's gap must be equal to : 

M = -(Amp - Ap)E = 2itr p (r p - v r )E 

where £ is the unperturbed gas density, to be identified with Elp74 above. The mass flow 
into the inner disc through the planet's gap corresponds to the angular momentum loss 
of the outer disc. It can be expressed as : 

M' = -j /j (16) 

where J Q is the total angular momentum of the outer disc and j is the specific angular 
momentum. Obviously, J Q = J ,un P — Jp, where J ,un P is the change of angular momentum 
of the outer disc in the unperturbed case and J p is the change of angular momentum of 
the planet. 

Thus, the mass flow deficit into the inner disc is : 

SM = M - M' = 2nr p (r p - v r )E + ( j GiUnp - j p )/j . (17) 

Now, remembering that J p = f p M p /2^/r^ and that j 0iU np — 27rr p i; r £j, the expression 
above becomes : 

«=f(^"l), (18) 

where \x = irr 2 Yi/M* and q = M p /M^ are the reduced masses of the inner disc and the 
planet. 

Some remarks on equation (ITS]) are in order. If the planet does not migrate, J p = 
and thus there is no mass flow deficit into the inner disc, in agreement with Fig. [HJ If 
q = 4/x there is also no mass deficit in the inner disc. This case corresponds to f p = v r , 
that is to a planet that migrates inward at the same speed of the unperturbed gas. We 
recall that this is the maximal migration speed of a planet in type II migration. So, if 
q < 4/i, equation (ITS]) does not apply. 

Finally, to compute the mass of the inner part of the disc (and therefore the cavity 
depth) one can follow the approach illustrated in Section [5] and multiply the result by : 

o/4u - 1 

1 - SM/M = 1 H -^— . (19) 

Q <^o,unp j 

4/i j p 
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Figure 14: Gap profile in a small mass evolving disc with 1Z = 10 5 and H/r = 0.03: 
the thin solid curve results from a numerical simulation while the bold line comes from 
our simple model. Thin dotted curve: the unperturbed profile £lp74(?", £„/2) from equa- 
tion ([2]). Bold dotted curve: the inner disc profile in our model, multiplied by the 
estimated relative mass deficit. 

From equation ( fl9j) . it also appears that if the planet mass is large with respect to the 
disc mass, then the mass flow deficit SM/Mq tends to J P /Jo, U np- If this planet opens a 
clean gap, J p ~ J ,unp an d r' p /v r <C 1, so that 5M « M . In summary, the two conditions 
q/fi ^> 1 and formation of a clean gap are required to deplete strongly the inner disc and 
form a deep cavity. 

To confirm this result, we performed a simulation with the same Jupiter mass planet 
initially placed on a circular orbit at r p = 5au, and the same initial disc profile as before, 
except that the gas density is divided by 100. The Reynolds number is taken equal to 
10 5 and the aspect ratio is decreased to 0.03 for the planet to open a cleaner gap than 
before. Fig. [TH displays the resulting disc profile at t = t u /2 (thin solid curve) : one 
can see that the planet did almost not migrate, but the inner disc has been significantly 
depleted with respect to the unperturbed profile (thin dotted curve) and to the profile 
given by equation (J3J) (bold solid curve). As we have seen above, for 'reasonable' values 
of disc height and viscosity a Jupiter mass planet does not open a very clean gap. So, 
in practice, an effective cavity opening planet would need to have 3-5 Jupiter masses, 
embedded in a sub- Jovian mass disc. 
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The migration rate of the planet at t — t u /2 in the simulation is J p = 2.88 x 10 -11 . The 
unperturbed density is 8.48 x 1CT 8 M Q .au~ 2 so that q/Afi = 36.84, and the unperturbed 
radial velocity of the gas is —1.03 x 1CT 5 . Following equation (fl~9l) . this gives a relative 
mass deficit in the inner disc : 

5M/M = 0.45 . 

Our model density profile, multiplied by 0.55 gives for the inner disc the bold dotted 
curve in Fig. [TH This agreement is quite satisfactory, although -admittedly- a better 
agreement with the real inner disc profile would be given by a coefficient of about 0.75. 



8 Gas cavities versus dust cavities 

Up to this point, we have discussed only the gas distribution of the disc. However, the 
observations reviewed in Section [2] constrain the absence of dust in the inner part of the 
disc, rather than the absence of gas. 

The dust dis tribution and the gas distrib ution are not necessarily the sam e. As recently 
pointed out by iRice et al. fbOQfil ) (see also IPaardekooper fc MellemalbOPfil ). the opening 



of a gap in the gas disc can act as a filter on the dust. Only the very small dust can 
flow with the gas through the gap and refill the inner part of the disc. Larger dust 
particles are repelled by the positive gas pressure gradient at the outer edge of the gap, 
and consequently they cannot pass through the planet's orbit. 

Our model provides all the ingredients to compute what is the maximal size of dust 
particles that can flow with the gas through the gap as a function of the various parameters 
of the problem (planet mass, V etc.), as outlined below. 

The gas orbital velocity v g differs from the Keplerian velocity vk due to pressure effects 

as: 

r r £ dr ' 

where P is the pressure and £ is the surface density. Using our equation of state (see 
Section 13.21) and assuming that sound speed c s does not change significantly along the 
gap's edge, the last term of the above equation can be rewritten as (l/£)d£/dr. Its 
maximum value can be obtained from equation ([3]) at r = r p + 2Rh- 

Once v g is computed, the outward drift speed of a dust particle of a given size and 
density can be computed following Weidenschilling ( 1977 ). The quantity Vr is actually a 



speed relative to the gas radial motion. Thus, if the inward radial velocity of the gas vi 9 ^ 
is larger in absolute value than Vr , the dust flows inwards, passing through the planet's 
gap. In the opposite case the dust is repelled by the pressure gradient and accumulates 
at the edge of the gap. The radial velocity of the gas at the base of the gap, where the 
pressure effect on the dust is maximal, can be computed as 

v (9) = ML (20) 

r 2vrr p £/(P) W 

where M' is given in equation (JTBl) . Assuming that the mass flux from the outer to the 
inner disc is not strongly perturbed by the presence of the planet (i.e. M' ~ Mo, which 
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is strictly true in the case of a non-migrating planet) one has : 



«?' = (21) 

where v r is given in equation Finally, because depends on the dust particle size, 
the equation \vj-\ = \vi 9 ^\ gives the maximal size of the particles that flow in the inner 
disc with the gas. 

However, we think that it is not evident that a real dust cavity is opened inside the 
orbit of the planet, even for dust sizes that cannot pass through the gap. The dust in 
the inner disc should accumulate at the place of the relative maximum of the gas density 
distribution, rather than follow the gas in its accretional motion on to the star. In fact, 
the positive density gradient at the inner edge of the disc is comparable to that at the 
outer edge of the gap opened by the planet (see Figs [T] and [7j) . So, if the dust is repelled 
by the gap edge, it should also be repelled by the inner disc edge. In this situation a ring 
of dust should form in the inner disc, its width depending on the local velocity dispersion 
due to the tu rbulence. Hence, the planet app ears opening a wide and deep gap in the dust 



distribution (IPaardekooper fc Mellemall2006l ). rather than a cavity. This, however, may 



still be consistent with the observations : if these rings of large dust grains are very narrow 
(as they may well be) the particles may not contribute significantly to the observed SED ; 
so, the SED observations would still suggest an inner cavity with small grains (producing 
a 10 fim feature) and a population extending to larger sizes in the outer disc. 

Possibilities to get rid of the larger dust grains in the inner disc could be (i) the 
collisional comminution into smaller particles that are coupled to the gas and therefore 
- as the gas - accrete on to the star, (ii) formation of larger particles or planetesimals, 
favoured by the high dust density in the ring and (iii) sublimation of the dust if the ring 
is located sufficiently close to the star. Obviously the situation is complicated. The study 
of the dust distribution requires the use of an adapted bi-fluid code, including irradiation 
effects. This should be the object of future work. 



9 Conclusion 

In this paper, we have performed numerical simulations of giant planets embedded in a 
gas disc of equivalent mass. The disc in the simulations is viscously evolving, reproducing 
(in absence of planetary perturbation) the evolution described in LP74. This is made 
possible by the use of a new numerical scheme that complements the usual 2D grid wit h 



a ID grid spanning over the assumed physical extension of the disc ( jCrida et al.l 120071 ). 
so that the global evolution of the disc is taken into account. This is a key point for the 
accurate simulation of both the cavity opening process and of planet migration. Indeed, 
the first phenomenon depends on the accretion of the inner disc on to the central star ; 
thus the simulation of the disc evolution down to its physical inner edge is necessary for 
reliable results. As for the planetary migration, proper type II migration corresponds to 
the case where the planet follows the disc evolution ; the global evolution of the disc is 
thus the key phenomenon that needs to be reproduced. 
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If the disc mass is not negligible with respect to the planet mass, we found that the 
planet does not modify substantially the disc profile with respect to the free viscous evo- 
lution ; it simply opens a more or less deep gap in the unperturbed profile. Consequently 
the opening of a cavity mainly depends on the position of the planetary orbit with respect 
to the inner edge of the disc R in {. Indeed, at this inner edge, the density is 0. It grows 
with r until the inner edge of the gap. Thus, if the gap is close to the inner edge of the 
disc, the density in the inner disc does not reach the maximum value and a cavity appears. 
As a result, it is much easier to open cavities is discs with large i?i n f than in discs with 
small i?i n f. The size of this inner radius may vary in astrophysical discs, depending on the 
process that governs the fall of the gas on to the central star. We stress that one should 
be aware of the primordial influence of this parameter on the inner disc evolution when 
computing numerical simulations. 

In the intermediate regime between type I and type II migration - where the planet 
opens a non-gas-proof gap- the migration may be stopped or reversed. In fact, if the 
gap is not clean, the gas in the gap exerts a positive viscous torque on the outer disc, 
so that the outer disc doesn't push the planet inward as efficiently as when the gap is 
clear. In addition, the corotation torque is positive and proportional to the density in the 
corotating zone. Thus, the total torque felt by the planet may be zero, or even positive, if 
the Reynolds number is low enough (Fig. [9]). We built a model based on simple, qualitative 
ideas, that leads to a simple expression of the total torque felt by a planet in an evolving 
disc (equation ( fl5i) ). This enabled us to reproduce the dependence of planet migration on 
the various parameters. Viscosity plays a major role : standard type II migration occurs 
only when the planet opens a very clean gap, that is at low viscosity ; if the viscosity 
increases, the gap becomes less depleted and less gas-proof, which leads to the decoupling 
of the migration with respect to the disc evolution. The role of the planet mass is very 
intuitive: the more massive is the planet, the deeper is its gap and thus the closer to 
standard type II migration is its behaviour. Last, the role of the radius of inner edge of 
the disc i?i n f is significant only when it is more than about half the orbital radius of the 
planet : in that case, the density gradient of the disc at the planet location increases with 
i?i n f, which enhances the corotation torque and makes the planet migrate outward. 

In conclusion, under some conditions on the disc parameters, type II migration may 
be avoided for a Jupiter mass planet at 5 au in an accreting disc, provided it does not 
open a very deep gap. This could explain why all the known giant planets (in the solar 
system and in extra-solar systems) are not hot Jupiters. Our results also explain why hot 
Jupiters had to stop migrating before falling on to their parent stars. They also explain 
why, among pairs of resonant planets, the outermost one is typically the most massive 
object. 
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